NASA Earthdata Cloud Clinic

Summary

Welcome to the NASA Earthdata Cloud Clinic!

We will focus on NASA Earthdata search & access.

We will use earthaccess for data search and direct cloud access, followed by xarray` for subsetting. Both are open source python libraries. We will also discover data using Earthdata Search

We will be accessing data directly from Amazon Web Services (AWS), specifically in the us-west-2 region, which is where all cloud-hosted NASA Earthdata reside. This shared compute environment (JupyterHub) is also running in the same location. We will then load the data into Python as an xarray dataset.

See the bottom of the notebook for additional resources, including several tutorials that served as a foundation for this clinic.

A note on earthaccess python library

In this example we will use the earthaccess library to search for data collections from NASA Earthdata. earthaccess is a Python library that simplifies data discovery and access to NASA Earth science data by providing an abstraction layer for NASA’s Common Metadata Repository (CMR) API Search API. The library makes searching for data more approachable by using a simpler notation instead of low level HTTP queries. earthaccess takes the trouble out of Earthdata Login authentication, makes search easier, and provides a stream-line way to download or stream search results into an xarray object for easy data access. It can be used on and off the cloud.

For more on earthaccess visit the earthaccess documentation site. Be aware that earthaccess is under active development, and your use and feedback help improve it!

A note on subsetting

In addition to directly accessing the files archived and distributed by each of the NASA DAACs, many datasets also support services that allow us to customize the data via subsetting, reformatting, reprojection/regridding, and file aggregation. What does subsetting mean? Here’s a generalized graphic of what we mean.

A note on jargon:

“direct cloud access” goes by several other names including “direct S3 access”, “direct access”, “direct in-region access”, “in-cloud data streaming”. And “subsetting” is also called “transformation”.

Learning Objectives

1. Utilize the earthaccess python library to search for data using spatial and temporal filters and explore search results
2. Stream data (i.e. perform in-region direct access of data) from an Amazon Simple Storage Service (S3) bucket where NASA Earthdata data is archived into our own cloud workspace, here in the Jupyter Lab/Notebook.
3. Extract variables and spatial slices from an xarray dataset

4. Plot data
   

5. Conceptualize data subsetting services provided by NASA Earthdata, including Harmony
6. Plot a polygon geojson file with a basemap using geoviews

---

Prerequisites

First we’ll import python packages and set our authentication that will be used for both of our access and subsetting methods.

This tutorial is meant to be run in the AWS cloud in the us-west-2 region. You’ll need to be aware that data in NASA’s Earthdata Cloud reside in Amazon Web Services (AWS) Simple Storage Service (S3) buckets. Access is provided via temporary credentials; this free access is limited to requests made within the US West (Oregon) (code: us-west-2) AWS region. While this compute location is required for direct S3 access, all data in Earthdata Cloud are still freely available via download.

Import Required Packages

# Suppress warnings
import warnings
warnings.simplefilter('ignore')
warnings.filterwarnings('ignore')

# Search & Access
import earthaccess 
from pprint import pprint
import xarray as xr
import hvplot.xarray #plot

# Harmony services
import geopandas as gpd
import geoviews as gv
gv.extension('bokeh', 'matplotlib')
from harmony import BBox, Client, Collection, Request, LinkType
import datetime as dt
import s3fs
%matplotlib inline

Authentication for NASA Earthdata

An Earthdata Login account is required to access data from the NASA Earthdata system. If you don’t already have one, visit https://urs.earthdata.nasa.gov to register and manage your Earthdata Login account. This account is free to create and only takes a moment to set up.

The first step is to get the correct authentication that will allow us to get cloud-hosted data from NASA. This is all done through Earthdata Login. We can use the earthaccess library here, where the login method also gets the correct AWS credentials.

The first time you run this it will ask you for your Earthdata Login username and password, and stores it in a .netrc file. After that, each time you authenticate with auth = earthaccess.login() it will log you in automatically.

# auth = earthaccess.login(strategy="interactive", persist=True)
auth = earthaccess.login()

1. earthaccess + xarray

earthaccess python library is an open-source library to simplify Earthdata Cloud search and access.

Search for data

There are multiple keywords we can use to discovery data from collections such as short_name, concept_id, and doi. The table below contains the short_name for some collections we are interested in for other exercises. Each of these can be used to search for data or information related to the collection we are interested in.

Shortname	Description	Example Temporal/Spatial parameters
SEA_SURFACE_HEIGHT_ALT_GRIDS_L4_2SATS_5DAY_6THDEG_V_JPL2205	enter description	temporal=(“2021-07-01”, “2021-09-30”)
MUR25-JPL-L4-GLOB-v04.2	MUR Sea Surface Temperature	temporal=(“2023-07-01”, “2023-07-03”), bounding_box=(-99, 18.19232, -78.85547, 31.23754)
EMITL2BMIN	EMIT L2B estimated mineral identification and band depths in a spatially raw, non-orthocorrected format	temporal=(“2023-07-01”, “2023-07-31”), bounding_box=(-99, 18.19232, -78.85547, 31.23754)
ECO_L2G_LSTE	ECOSTRESS Gridded Land Surface Temperature and Emissivity Instantaneous Level 2 Global 70 m	temporal=(“2023-07-01”, “2023-07-03”), bounding_box=(-99, 18.19232, -78.85547, 31.23754)

But wait…You may be asking “how can we find the short_name for collections not in the table above?”.

–> Let’s take a quick detour and head to Earthdata Search GUI to gather more information about our dataset of interest. The dataset “short name” can be found by clicking on the Info button on our collection search result, and we can paste that into a python variable.

(Side Note: Both earthaccess python libray and the Earthdata Search (GUI) leverage the Common Metadata Repository (CMR) API to search for collections and granules.)

Here we use the search_data function of earthaccess to query based on the short_name of interest, as well as other paramters such as temporal range:

data_name = "SEA_SURFACE_HEIGHT_ALT_GRIDS_L4_2SATS_5DAY_6THDEG_V_JPL2205"

results = earthaccess.search_data(
    short_name=data_name,
    cloud_hosted=True,
    temporal=("2021-07-01", "2021-09-30"),
)

Granules found: 18

According to PO.DAAC’s dataset landing page, gridded Sea Surface Height Anomalies (SSHA) above a mean sea surface are provided. The data are produced on a 1/6th degree grid every 5 days.

We can discover more information about the matching files:

pprint(results[0])

Collection: {'Version': '2205', 'ShortName': 'SEA_SURFACE_HEIGHT_ALT_GRIDS_L4_2SATS_5DAY_6THDEG_V_JPL2205'}
Spatial coverage: {'HorizontalSpatialDomain': {'Geometry': {'BoundingRectangles': [{'WestBoundingCoordinate': 0.083, 'SouthBoundingCoordinate': -79.917, 'EastBoundingCoordinate': 180, 'NorthBoundingCoordinate': 79.917}, {'WestBoundingCoordinate': -180, 'SouthBoundingCoordinate': -79.917, 'EastBoundingCoordinate': -0.083, 'NorthBoundingCoordinate': 79.917}]}}}
Temporal coverage: {'RangeDateTime': {'EndingDateTime': '2021-07-05T00:00:00.000Z', 'BeginningDateTime': '2021-07-05T00:00:00.000Z'}}
Size(MB): 9.307239532470703
Data: ['https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/SEA_SURFACE_HEIGHT_ALT_GRIDS_L4_2SATS_5DAY_6THDEG_V_JPL2205/ssh_grids_v2205_2021070512.nc']

Access data

Our code will work the same way if we are running it “in-region”, within our shared cloud environment, or locally from our laptop.

Since we are working in the AWS us-west-2 region, we can stream data directly to xarray. We are using the open_mfdataset() (multi-file) method, which is required when using earthaccess.

(Tips: To open a single file, if troubleshooting for example: ds = xr.open_dataset(earthaccess.open(results)[0]))

ds = xr.open_mfdataset(earthaccess.open(results))
ds

Opening 18 granules, approx size: 0.16 GB
using endpoint: https://archive.podaac.earthdata.nasa.gov/s3credentials

<xarray.Dataset> Size: 299MB
Dimensions:      (Time: 18, Longitude: 2160, nv: 2, Latitude: 960)
Coordinates:
  * Longitude    (Longitude) float32 9kB 0.08333 0.25 0.4167 ... 359.8 359.9
  * Latitude     (Latitude) float32 4kB -79.92 -79.75 -79.58 ... 79.75 79.92
  * Time         (Time) datetime64[ns] 144B 2021-07-05T12:00:00 ... 2021-09-2...
Dimensions without coordinates: nv
Data variables:
    Lon_bounds   (Time, Longitude, nv) float32 311kB dask.array<chunksize=(1, 2160, 2), meta=np.ndarray>
    Lat_bounds   (Time, Latitude, nv) float32 138kB dask.array<chunksize=(1, 960, 2), meta=np.ndarray>
    Time_bounds  (Time, nv) datetime64[ns] 288B dask.array<chunksize=(1, 2), meta=np.ndarray>
    SLA          (Time, Latitude, Longitude) float32 149MB dask.array<chunksize=(1, 960, 2160), meta=np.ndarray>
    SLA_ERR      (Time, Latitude, Longitude) float32 149MB dask.array<chunksize=(1, 960, 2160), meta=np.ndarray>
Attributes: (12/21)
    Conventions:            CF-1.6
    ncei_template_version:  NCEI_NetCDF_Grid_Template_v2.0
    Institution:            Jet Propulsion Laboratory
    geospatial_lat_min:     -79.916664
    geospatial_lat_max:     79.916664
    geospatial_lon_min:     0.083333336
    ...                     ...
    version_number:         2205
    Data_Pnts_Each_Sat:     {"16": 743215, "1007": 674076}
    source_version:         commit 58c7da13c0c0069ae940c33a82bf1544b7d991bf
    SLA_Global_MEAN:        0.06428374482174487
    SLA_Global_STD:         0.0905195660534004
    latency:                final

xarray.Dataset

• Dimensions:
  • Time: 18
  • Longitude: 2160
  • nv: 2
  • Latitude: 960
• Coordinates: (3)
  • Longitude
    (Longitude)
    float32
    0.08333 0.25 0.4167 ... 359.8 359.9

    standard_name :

    longitude

    units :

    degrees_east

    point_spacing :

    even

    long_name :

    longitude

    axis :

    X

    bounds :

    Lon_bounds

    array([8.333334e-02, 2.500000e-01, 4.166667e-01, ..., 3.595833e+02,
           3.597500e+02, 3.599167e+02], dtype=float32)

  • Latitude
    (Latitude)
    float32
    -79.92 -79.75 ... 79.75 79.92

    standard_name :

    latitude

    units :

    degrees_north

    point_spacing :

    even

    long_name :

    latitude

    axis :

    Y

    bounds :

    Lat_bounds

    array([-79.916664, -79.75    , -79.583336, ...,  79.583336,  79.75    ,
            79.916664], dtype=float32)

  • Time
    (Time)
    datetime64[ns]
    2021-07-05T12:00:00 ... 2021-09-...

    standard_name :

    time

    long_name :

    Time

    bounds :

    Time_bounds

    axis :

    T

    array(['2021-07-05T12:00:00.000000000', '2021-07-10T12:00:00.000000000',
           '2021-07-15T12:00:00.000000000', '2021-07-20T12:00:00.000000000',
           '2021-07-25T12:00:00.000000000', '2021-07-30T12:00:00.000000000',
           '2021-08-04T12:00:00.000000000', '2021-08-09T12:00:00.000000000',
           '2021-08-14T12:00:00.000000000', '2021-08-19T12:00:00.000000000',
           '2021-08-24T12:00:00.000000000', '2021-08-29T12:00:00.000000000',
           '2021-09-03T12:00:00.000000000', '2021-09-08T12:00:00.000000000',
           '2021-09-13T12:00:00.000000000', '2021-09-18T12:00:00.000000000',
           '2021-09-23T12:00:00.000000000', '2021-09-28T12:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (5)
  • Lon_bounds
    (Time, Longitude, nv)
    float32
    dask.array<chunksize=(1, 2160, 2), meta=np.ndarray>

    units :

    degrees_east

    comment :

    longitude values at the west and east bounds of each pixel.

    Array Chunk Bytes 303.75 kiB 16.88 kiB Shape (18, 2160, 2) (1, 2160, 2) Dask graph 18 chunks in 55 graph layers Data type float32 numpy.ndarray

    Lat_bounds

    (Time, Latitude, nv)

    float32

    dask.array<chunksize=(1, 960, 2), meta=np.ndarray>

    

    units :

    degrees_north

    comment :

    latitude values at the north and south bounds of each pixel.

    Array Chunk Bytes 135.00 kiB 7.50 kiB Shape (18, 960, 2) (1, 960, 2) Dask graph 18 chunks in 55 graph layers Data type float32 numpy.ndarray

    Time_bounds

    (Time, nv)

    datetime64[ns]

    dask.array<chunksize=(1, 2), meta=np.ndarray>

    

    comment :

    Time bounds for each time value, same value as time variable. The time variable is defined on points instead of on bounding boxes.

    Array Chunk Bytes 288 B 16 B Shape (18, 2) (1, 2) Dask graph 18 chunks in 37 graph layers Data type datetime64[ns] numpy.ndarray

    SLA

    (Time, Latitude, Longitude)

    float32

    dask.array<chunksize=(1, 960, 2160), meta=np.ndarray>

    

    units :

    m

    long_name :

    Sea Level Anomaly Estimate

    standard_name :

    sea_surface_height_above_sea_level

    alias :

    sea_surface_height_above_sea_level

    Array Chunk Bytes 142.38 MiB 7.91 MiB Shape (18, 960, 2160) (1, 960, 2160) Dask graph 18 chunks in 37 graph layers Data type float32 numpy.ndarray

    SLA_ERR

    (Time, Latitude, Longitude)

    float32

    dask.array<chunksize=(1, 960, 2160), meta=np.ndarray>

    

    units :

    m

    long_name :

    Sea Level Anomaly Error Estimate

    Array Chunk Bytes 142.38 MiB 7.91 MiB Shape (18, 960, 2160) (1, 960, 2160) Dask graph 18 chunks in 37 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • Longitude
    PandasIndex

    PandasIndex(Index([0.0833333358168602,               0.25, 0.4166666567325592,
           0.5833333134651184,               0.75, 0.9166666865348816,
           1.0833333730697632,               1.25, 1.4166666269302368,
           1.5833333730697632,
           ...
            358.4166564941406,  358.5833435058594,             358.75,
            358.9166564941406,  359.0833435058594,             359.25,
            359.4166564941406,  359.5833435058594,             359.75,
            359.9166564941406],
          dtype='float32', name='Longitude', length=2160))

  • Latitude
    PandasIndex

    PandasIndex(Index([-79.91666412353516,             -79.75, -79.58333587646484,
           -79.41666412353516,             -79.25, -79.08333587646484,
           -78.91666412353516,             -78.75, -78.58333587646484,
           -78.41666412353516,
           ...
            78.41666412353516,  78.58333587646484,              78.75,
            78.91666412353516,  79.08333587646484,              79.25,
            79.41666412353516,  79.58333587646484,              79.75,
            79.91666412353516],
          dtype='float32', name='Latitude', length=960))

  • Time
    PandasIndex

    PandasIndex(DatetimeIndex(['2021-07-05 12:00:00', '2021-07-10 12:00:00',
                   '2021-07-15 12:00:00', '2021-07-20 12:00:00',
                   '2021-07-25 12:00:00', '2021-07-30 12:00:00',
                   '2021-08-04 12:00:00', '2021-08-09 12:00:00',
                   '2021-08-14 12:00:00', '2021-08-19 12:00:00',
                   '2021-08-24 12:00:00', '2021-08-29 12:00:00',
                   '2021-09-03 12:00:00', '2021-09-08 12:00:00',
                   '2021-09-13 12:00:00', '2021-09-18 12:00:00',
                   '2021-09-23 12:00:00', '2021-09-28 12:00:00'],
                  dtype='datetime64[ns]', name='Time', freq=None))

• Attributes: (21)

  Conventions :

  CF-1.6

  ncei_template_version :

  NCEI_NetCDF_Grid_Template_v2.0

  Institution :

  Jet Propulsion Laboratory

  geospatial_lat_min :

  -79.916664

  geospatial_lat_max :

  79.916664

  geospatial_lon_min :

  0.083333336

  geospatial_lon_max :

  359.91666

  time_coverage_start :

  2021-07-05

  time_coverage_end :

  2021-07-05

  date_created :

  2022-10-30T20:58:16.297288

  title :

  Sea Level Anomaly Estimate based on Altimeter Data, final product (replaced interim version).

  short_name :

  SEA_SURFACE_HEIGHT_ALT_GRIDS_L4_2SATS_5DAY_6THDEG_V_JPL2205

  long_name :

  MEaSUREs Gridded Sea Surface Height Anomalies Version 2205

  summary :

  Sea level anomaly grids from altimeter data using Kriging interpolation, which gives best linear prediction based upon prior knowledge of covariance.

  DOI :

  10.5067/SLREF-CDRV3

  version_number :

  2205

  Data_Pnts_Each_Sat :

  {"16": 743215, "1007": 674076}

  source_version :

  commit 58c7da13c0c0069ae940c33a82bf1544b7d991bf

  SLA_Global_MEAN :

  0.06428374482174487

  SLA_Global_STD :

  0.0905195660534004

  latency :

  final

Plot the data

Let’s make a quick interactive plot of the data using an open source tool call hvplot. Because our data is 3d and has a time component, we can also preview the data over time, using the slider of hvplot.

ds.SLA.hvplot.image(x='Longitude', y='Latitude', cmap='RdBu', clim=(-2,2), title="Sea Level Anomaly Estimate (m)")

Subset data

We can create a subsetted xarray dataset by extracting the SLA variable and slicing the dataset by a smaller area of interest near the state of Baja, Mexico.

ds_subset = ds['SLA'].sel(Latitude=slice(15.8, 35.9), Longitude=slice(234.5,260.5)) 
ds_subset

<xarray.DataArray 'SLA' (Time: 18, Latitude: 120, Longitude: 156)> Size: 1MB
dask.array<getitem, shape=(18, 120, 156), dtype=float32, chunksize=(1, 120, 156), chunktype=numpy.ndarray>
Coordinates:
  * Longitude  (Longitude) float32 624B 234.6 234.8 234.9 ... 260.1 260.2 260.4
  * Latitude   (Latitude) float32 480B 15.92 16.08 16.25 ... 35.42 35.58 35.75
  * Time       (Time) datetime64[ns] 144B 2021-07-05T12:00:00 ... 2021-09-28T...
Attributes:
    units:          m
    long_name:      Sea Level Anomaly Estimate
    standard_name:  sea_surface_height_above_sea_level
    alias:          sea_surface_height_above_sea_level

xarray.DataArray

'SLA'

• Time: 18
• Latitude: 120
• Longitude: 156

• dask.array<chunksize=(1, 120, 156), meta=np.ndarray>

  Array Chunk Bytes 1.29 MiB 73.12 kiB Shape (18, 120, 156) (1, 120, 156) Dask graph 18 chunks in 38 graph layers Data type float32 numpy.ndarray

• Coordinates: (3)
  • Longitude
    (Longitude)
    float32
    234.6 234.8 234.9 ... 260.2 260.4

    standard_name :

    longitude

    units :

    degrees_east

    point_spacing :

    even

    long_name :

    longitude

    axis :

    X

    bounds :

    Lon_bounds

    array([234.58333, 234.75   , 234.91667, 235.08333, 235.25   , 235.41667,
           235.58333, 235.75   , 235.91667, 236.08333, 236.25   , 236.41667,
           236.58333, 236.75   , 236.91667, 237.08333, 237.25   , 237.41667,
           237.58333, 237.75   , 237.91667, 238.08333, 238.25   , 238.41667,
           238.58333, 238.75   , 238.91667, 239.08333, 239.25   , 239.41667,
           239.58333, 239.75   , 239.91667, 240.08333, 240.25   , 240.41667,
           240.58333, 240.75   , 240.91667, 241.08333, 241.25   , 241.41667,
           241.58333, 241.75   , 241.91667, 242.08333, 242.25   , 242.41667,
           242.58333, 242.75   , 242.91667, 243.08333, 243.25   , 243.41667,
           243.58333, 243.75   , 243.91667, 244.08333, 244.25   , 244.41667,
           244.58333, 244.75   , 244.91667, 245.08333, 245.25   , 245.41667,
           245.58333, 245.75   , 245.91667, 246.08333, 246.25   , 246.41667,
           246.58333, 246.75   , 246.91667, 247.08333, 247.25   , 247.41667,
           247.58333, 247.75   , 247.91667, 248.08333, 248.25   , 248.41667,
           248.58333, 248.75   , 248.91667, 249.08333, 249.25   , 249.41667,
           249.58333, 249.75   , 249.91667, 250.08333, 250.25   , 250.41667,
           250.58333, 250.75   , 250.91667, 251.08333, 251.25   , 251.41667,
           251.58333, 251.75   , 251.91667, 252.08333, 252.25   , 252.41667,
           252.58333, 252.75   , 252.91667, 253.08333, 253.25   , 253.41667,
           253.58333, 253.75   , 253.91667, 254.08333, 254.25   , 254.41667,
           254.58333, 254.75   , 254.91667, 255.08333, 255.25   , 255.41667,
           255.58333, 255.75   , 255.91667, 256.08334, 256.25   , 256.41666,
           256.58334, 256.75   , 256.91666, 257.08334, 257.25   , 257.41666,
           257.58334, 257.75   , 257.91666, 258.08334, 258.25   , 258.41666,
           258.58334, 258.75   , 258.91666, 259.08334, 259.25   , 259.41666,
           259.58334, 259.75   , 259.91666, 260.08334, 260.25   , 260.41666],
          dtype=float32)

  • Latitude
    (Latitude)
    float32
    15.92 16.08 16.25 ... 35.58 35.75

    standard_name :

    latitude

    units :

    degrees_north

    point_spacing :

    even

    long_name :

    latitude

    axis :

    Y

    bounds :

    Lat_bounds

    array([15.916667, 16.083334, 16.25    , 16.416666, 16.583334, 16.75    ,
           16.916666, 17.083334, 17.25    , 17.416666, 17.583334, 17.75    ,
           17.916666, 18.083334, 18.25    , 18.416666, 18.583334, 18.75    ,
           18.916666, 19.083334, 19.25    , 19.416666, 19.583334, 19.75    ,
           19.916666, 20.083334, 20.25    , 20.416666, 20.583334, 20.75    ,
           20.916666, 21.083334, 21.25    , 21.416666, 21.583334, 21.75    ,
           21.916666, 22.083334, 22.25    , 22.416666, 22.583334, 22.75    ,
           22.916666, 23.083334, 23.25    , 23.416666, 23.583334, 23.75    ,
           23.916666, 24.083334, 24.25    , 24.416666, 24.583334, 24.75    ,
           24.916666, 25.083334, 25.25    , 25.416666, 25.583334, 25.75    ,
           25.916666, 26.083334, 26.25    , 26.416666, 26.583334, 26.75    ,
           26.916666, 27.083334, 27.25    , 27.416666, 27.583334, 27.75    ,
           27.916666, 28.083334, 28.25    , 28.416666, 28.583334, 28.75    ,
           28.916666, 29.083334, 29.25    , 29.416666, 29.583334, 29.75    ,
           29.916666, 30.083334, 30.25    , 30.416666, 30.583334, 30.75    ,
           30.916666, 31.083334, 31.25    , 31.416666, 31.583334, 31.75    ,
           31.916666, 32.083332, 32.25    , 32.416668, 32.583332, 32.75    ,
           32.916668, 33.083332, 33.25    , 33.416668, 33.583332, 33.75    ,
           33.916668, 34.083332, 34.25    , 34.416668, 34.583332, 34.75    ,
           34.916668, 35.083332, 35.25    , 35.416668, 35.583332, 35.75    ],
          dtype=float32)

  • Time
    (Time)
    datetime64[ns]
    2021-07-05T12:00:00 ... 2021-09-...

    standard_name :

    time

    long_name :

    Time

    bounds :

    Time_bounds

    axis :

    T

    array(['2021-07-05T12:00:00.000000000', '2021-07-10T12:00:00.000000000',
           '2021-07-15T12:00:00.000000000', '2021-07-20T12:00:00.000000000',
           '2021-07-25T12:00:00.000000000', '2021-07-30T12:00:00.000000000',
           '2021-08-04T12:00:00.000000000', '2021-08-09T12:00:00.000000000',
           '2021-08-14T12:00:00.000000000', '2021-08-19T12:00:00.000000000',
           '2021-08-24T12:00:00.000000000', '2021-08-29T12:00:00.000000000',
           '2021-09-03T12:00:00.000000000', '2021-09-08T12:00:00.000000000',
           '2021-09-13T12:00:00.000000000', '2021-09-18T12:00:00.000000000',
           '2021-09-23T12:00:00.000000000', '2021-09-28T12:00:00.000000000'],
          dtype='datetime64[ns]')

• Indexes: (3)
  • Longitude
    PandasIndex

    PandasIndex(Index([234.5833282470703,            234.75, 234.9166717529297,
           235.0833282470703,            235.25, 235.4166717529297,
           235.5833282470703,            235.75, 235.9166717529297,
           236.0833282470703,
           ...
           258.9166564941406, 259.0833435058594,            259.25,
           259.4166564941406, 259.5833435058594,            259.75,
           259.9166564941406, 260.0833435058594,            260.25,
           260.4166564941406],
          dtype='float32', name='Longitude', length=156))

  • Latitude
    PandasIndex

    PandasIndex(Index([15.916666984558105,  16.08333396911621,              16.25,
            16.41666603088379,  16.58333396911621,              16.75,
            16.91666603088379,  17.08333396911621,              17.25,
            17.41666603088379,
           ...
                        34.25,  34.41666793823242,  34.58333206176758,
                        34.75,  34.91666793823242,  35.08333206176758,
                        35.25,  35.41666793823242,  35.58333206176758,
                        35.75],
          dtype='float32', name='Latitude', length=120))

  • Time
    PandasIndex

    PandasIndex(DatetimeIndex(['2021-07-05 12:00:00', '2021-07-10 12:00:00',
                   '2021-07-15 12:00:00', '2021-07-20 12:00:00',
                   '2021-07-25 12:00:00', '2021-07-30 12:00:00',
                   '2021-08-04 12:00:00', '2021-08-09 12:00:00',
                   '2021-08-14 12:00:00', '2021-08-19 12:00:00',
                   '2021-08-24 12:00:00', '2021-08-29 12:00:00',
                   '2021-09-03 12:00:00', '2021-09-08 12:00:00',
                   '2021-09-13 12:00:00', '2021-09-18 12:00:00',
                   '2021-09-23 12:00:00', '2021-09-28 12:00:00'],
                  dtype='datetime64[ns]', name='Time', freq=None))

• Attributes: (4)

  units :

  m

  long_name :

  Sea Level Anomaly Estimate

  standard_name :

  sea_surface_height_above_sea_level

  alias :

  sea_surface_height_above_sea_level

ds_subset.hvplot.image(x='Longitude', y='Latitude', cmap='RdBu', clim=(-0.5, 0.5), title="Sea Level Anomaly Estimate (m)")

---

Exercise

Let’s explore other datasets from the table above

MUR SST

data_name = "MUR25-JPL-L4-GLOB-v04.2"

sst_results = earthaccess.search_data(
    short_name=data_name,
    cloud_hosted=True,
    temporal=("2023-07-01", "2023-07-03"),  
    #bounding_box(lower_left_lon, lower_left_lat, upper_right_lon, upper_right_lat) 
    bounding_box=(-99, 18.19232, -78.85547, 31.23754), 
)

pprint(sst_results[0])

ds = xr.open_mfdataset(earthaccess.open(sst_results))
ds

ds.sst_anomaly.hvplot.image(x='lon', y='lat', cmap='RdBu', clim=(-10,10), title="Sea Surfate Temperature Anomaly (K)")

ECOSTRESS

data_name = "ECO_L2G_LSTE"

eco_results = earthaccess.search_data(
    short_name=data_name,
    cloud_hosted=True,
    temporal=("2023-07-01", "2023-07-03"),   # ECOSTRESS
    #bounding_box(lower_left_lon, lower_left_lat, upper_right_lon, upper_right_lat) 
    bounding_box=(-99, 18.19232, -78.85547, 31.23754),  # ECOSTRESS
)

pprint(results[0])

Download one file to do h5dump to see contents and groups, so we can open the right group in memory from the hdf5 file:

download:
download_files = earthaccess.download(results[0], "/home/jovyan/earthdata-cloud-cookbook/examples/data/")

h5dump:
h5dump -n 1 --contents=1 filename.h5 > dataname_h5dump.txt

file_handlers = earthaccess.open(eco_results)
file_handlers

DOES NOT WORK RIGHT NOW:

# Open the h5 file at a specific area using xarray, in memory
lst = xr.open_dataset(file_handlers[0], group='/HDFEOS/GRIDS/ECO_L2G_LSTE_70m/Data Fields/')
lst

---

EMIT

data_name = "EMITL2BMIN"

emit_results = earthaccess.search_data(
    short_name=data_name,
    cloud_hosted=True,
    temporal=("2023-07-01", "2023-07-31"),   
    #bounding_box(lower_left_lon, lower_left_lat, upper_right_lon, upper_right_lat) 
    bounding_box=(-99, 18.19232, -78.85547, 31.23754),
)

pprint(emit_results[0])

# ERROR when trying to open multiple EMIT granules: something about it can't concatenate 
ds = xr.open_mfdataset(earthaccess.open(emit_results))
ds

# Open one EMIT granule/file
ds = xr.open_dataset(earthaccess.open(emit_results)[0])
ds

Any suggestions for what and how to do a quick plot of an EMIT variable?

ds.group_1_mineral_id.hvplot.image(x='downtrack', y='crosstrack')

---

---

2. Harmony-py + xarray

Bonus cloud clinic topic that we won’t cover live together: subsetting with Earthdata Harmony API services

We can also use the Harmony-py python library for direct cloud access & data subsetting: * Requesting a subset of data from the GHRSST Level 4 MUR Global Foundation Sea Surface Temperature Analysis (v4.1) dataset using a vector-based geospatial file. * This dataset can also be viewed in Earthdata Search.

harmony-py + xarray

We will be accessing a subsetted NetCDF-4 file using Transformation Services in the Cloud.

Using Earthdata Search, we can find datasets that support these services using the “Customizable” filter.

We will find, request, and open customized data using Harmony, below.

Define area of interest

First, use geopandas to read in a geojson file containing a polygon feature over the Gulf of Mexico. The geojson file is found in the ~/data directory.

geojson_path = './data/gulf.json'
gdf = gpd.read_file(geojson_path) #Return a GeoDataFrame object

We can plot the polygon using the geoviews package that we imported as gv with ‘bokeh’ and ‘matplotlib’ extensions. The following has reasonable width, height, color, and line widths to view our polygon when it is overlayed on a base tile map. We can view this with the * operator.

base = gv.tile_sources.EsriImagery.opts(width=650, height=500)
ocean_map = gv.Polygons(gdf).opts(line_color='yellow', line_width=5, color=None)
base * ocean_map

Create a subset request

Here we’ll create a request for a spatial subset of data.

First, we need to create a Harmony Client, which is what we will interact with to submit and inspect a data request to Harmony, as well as to retrieve results.

When creating the Client, we need to provide Earthdata Login credentials. This basic line below assumes that we have a .netrc available. See the Earthdata Cloud Cookbook appendix for more information on Earthdata Login and netrc setup.

harmony_client = Client()

See the harmony-py documentation for details on how to construct your request.

sst_short_name="MUR-JPL-L4-GLOB-v4.1"

request = Request(
    collection=Collection(id=sst_short_name),
    shape=geojson_path,
    temporal={
    'start': dt.datetime(2021, 8, 1, 1),
    'stop': dt.datetime(2021, 8, 1, 2)   
    },
)

Submit a subset request

Now that the request is created, we can now submit it to Harmony using the Harmony Client object. A job id is returned, which is a unique identifier that represents the submitted request.

job_id = harmony_client.submit(request)
job_id

Check request status

Depending on the size of the request, it may be helpful to wait until the request has completed processing before the remainder of the code is executed. The wait_for_processing() method will block subsequent lines of code while optionally showing a progress bar.

harmony_client.wait_for_processing(job_id, show_progress=True)

We can view Harmony job response and output URLs Once the data request has finished processing, we can view details on the job that was submitted to Harmony, including the API call to Harmony, and informational messages on the request if available.

result_json() calls wait_for_processing() and returns the complete job in JSON format once processing is complete.

data = harmony_client.result_json(job_id)
pprint(data)

Access data

Just like above, the subsetted outputs produced by Harmony can be accessed directly from the cloud.

First, we’ll retrieve list of output URLs.

The result_urls() method calls wait_for_processing() and returns a list of the processed data URLs once processing is complete. You may optionally show the progress bar as shown below.

results = harmony_client.result_urls(job_id, link_type=LinkType.s3)
urls = list(results)
url = urls[0]
print(url)

Using aws_credentials you can retrieve the credentials needed to access the Harmony s3 staging bucket and its contents.

creds = harmony_client.aws_credentials()

Open staged files with s3fs and xarray

We use the AWS s3fs package to create a file system that can then be read by xarray:

s3_fs = s3fs.S3FileSystem(
    key=creds['aws_access_key_id'],
    secret=creds['aws_secret_access_key'],
    token=creds['aws_session_token'],
    client_kwargs={'region_name':'us-west-2'},
)

Now that we have our s3 file system set, including our declared credentials, we’ll use that to open the url, and read in the file through xarray. This extra step is needed because xarray cannot open the S3 location directly. Instead, the S3 file object is passed to xarray, in order to then open the dataset.

f = s3_fs.open(url, mode='rb')
ds = xr.open_dataset(f)
ds

As before, we use the xarray built in plotting function to create a simple plot along the x and y dimensions of the dataset. We can see that the data are subsetted to our polygon:

ds.analysed_sst.plot();

---

Additional Resources

Tutorials

This clinic was based off of several notebook tutorials including those presented during past workshop events, along with other materials co-created by the NASA Openscapes mentors: * 2021 Earthdata Cloud Hackathon * 2021 AGU Workshop * Accessing and working with ICESat-2 data in the cloud * Analyzing Sea Level Rise Using Earth Data in the Cloud

Cloud services

The examples used in the clinic provide an abbreviated and simplified workflow to explore access and subsetting options available through the Earthdata Cloud. There are several other options that can be used to interact with data in the Earthdata Cloud including:

• OPeNDAP
  • Hyrax provides direct access to subsetting of NASA data using Python or your favorite analysis tool
  • Tutorial highlighting OPeNDAP usage: https://nasa-openscapes.github.io/earthdata-cloud-cookbook/how-tos/working-locally/Earthdata_Cloud__Data_Access_OPeNDAP_Example.html
• Zarr-EOSDIS-Store
  • The zarr-eosdis-store library allows NASA EOSDIS Collections to be accessed efficiently by the Zarr Python library, provided they have a sidecar DMR++ metadata file generated.
  • Tutorial highlighting this library’s usage: https://nasa-openscapes.github.io/2021-Cloud-Hackathon/tutorials/09_Zarr_Access.html

Support

• Earthdata Forum
  • User Services and community support for all things NASA Earthdata, including Earthdata Cloud
• Earthdata Webinar series
  • Webinars from DAACs and other groups across EOSDIS including guidance on working with Earthdata Cloud
  • See the Earthdata YouTube channel for more videos